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ABSTRACT 


Cavities  generate  delayed  radar  returns.  But  because  the  standard  scattering  models  are 
overwhelmed  by  the  complexity  of  the  problem,  no  predictive  physical  model  of  cavity  scattering  has 
been  available  to  describe  what  is  observed. 

The  delayed  return  is  caused  by  the  slow  disgorgement  of  the  radar  energy  captured  by  the  cavity. 
For  some  cavities,  and  especially  for  the  liquid-fueled  rocket  engine,  the  length  of  the  extended  return  is 
many  times  the  maximum  dimension  of  the  cavity.  Since  this  implies  that  the  radar  energy  undergoes  a 
great  many  internal  reflections  before  it  reemerges,  it  suggests  treating  the  captured  energy  within  the 
cavity  as  a  well-mixed  thermal  radiation  field. 

This  suggestion  leads  to  a  simple  scattering  model  for  the  well-mixed  lossless  thermal  cavity  and  to 
a  straightforward  extension  to  accommodate  the  nozzle  on  the  engine.  The  model  gives,  independently  of 
both  polarization  and  frequency,  the  formula  G  =  2AF(Gi)F(6s )cos#(  cos#v  for  the  narrowband  bistatic 
RCS  of  the  rocket-engine’s  extended  signature,  w  here  0.  and  6?  are,  respectively,  the  angles  between  the 
engine  axis  and  the  lines  of  sight  to  the  transmitting  and  receiving  radars,  A  is  the  area  of  the  nozzle 
throat,  and  F(0)  is  the  nozzle  gain  factor.  For  a  cavity  with  no  nozzle,  F{6)  =  1 

The  concept  also  produces  a  somewhat  longer  formula  defining  the  wideband  extended  radar 
signature.  It  predicts,  independently  of  incidence  or  scattering  angles,  polarization  or  frequency,  a  decay 
rate  of  2.17/L*  dB/m.  ( L *  is  the  characteristic  chamber  length  defined  in  rocket  engineering  as  the  ratio  of 
the  chamber  volume  to  the  throat  area.) 

A  second  component  of  the  scattered  energy,  the  nonextended  return,  is  the  energy  reflected  by  the 
nozzle  without  entering  the  combustion  chamber.  Added  noncoherently  to  the  energy  disgorged  by  the 
combustion  chamber,  it  gives  the  total  signature. 

The  report  presents  the  derivation  of  these  formulas  and  then,  by  way  of  illustration,  applies  the 
model  to  a  hypothetical  rocket  engine.  It  also  includes  predictions  of  the  decay  rates  for  a  number  of 
rocket  engines,  both  foreign  and  domestic.  Remarkably,  the  decay  rates  of  six  different  high-thrust  rocket 
engines  all  lie  within  a  narrow  2-to-l  range,  whereas  their  thrusts  cover  a  180-to-l  range. 

The  energy  capture  and  release  processes  of  the  cavity  are  so  simple  to  model  that  the  overall  model 
accuracy  would  seem  to  depend  wholly  on  the  assumption  that  the  cavity  energy  is  well  mixed  a 
reasonably  safe  assumption  for  a  rocket  engine,  as  well  as  for  other  highly  cluttered  cavities.  However, 
the  report  includes  no  comparison  with  measured  data. 


This  page  intentionally  left  blank. 


ACKNOWLEDGMENTS 


I  would  like  to  acknowledge  the  encouragement  and  technical  stimulation  1  received  during  this 
work  from  Douglas  Koltenuk,  Ryan  Suess,  and  Allison  Duncan. 


v 


This  page  intentionally  left  blank. 


TABLE  OF  CONTENTS 


Abstract  in 

Acknowledgments  v 

List  of  Illustrations  ix 

1 .  CAVITY  SIGNATURE:  A  THERMAL  MODEL  1 

1.1  Introduction  1 

1.2  Basic  Thermal  Cavity  2 

1 .3  Nozzle  Gam  Factor  5 

1 .4  On-Axis  Gain  Factor  for  a  Conical  Nozzle  7 

1 .5  Discussion  8 

1.6  Conclusions  10 

APPENDIX  1  1 

A.  1  Basic  Well-Mixed  Thermal  Cavity  1 1 

A. 2  Wideband  Signature  1 1 

A. 3  Nozzle  Gain-Factor  Reciprocity  13 

A. 4  RCS  Lobe-Width  Estimate  13 

A. 5  On-Axis  Gam  Factor  for  a  Conical  Nozzle  14 

A. 6  Ray-Trace  MATLAB  Function  16 


vii 


This  page  intentionally  left  blank. 


LIST  OF  ILLUSTRATIONS 


Figure  Page 

No. 

1  The  basic  cavity,  showing  the  aperture,  the  single  incidence  direction,  and  the  2 

energy  discharging  along  multiple  scattering  directions. 

2  The  wideband  on-axis  backscattered  cavity  signature  for  the  basic  well-mixed  4 

thermal  cavity.  The  unrealistic  abrupt  start  of  the  original  formula  is  suppressed 

in  the  modified  formula.  The  y- axis  units  are  square  meters  of  radar  cross 
section  per  meter  of  range,  expressed  in  dB.  (The  effective  cavity  radius  is  the 
radius  of  a  sphere  having  the  same  volume  as  the  cavity.) 

3  Four  images  showing  ray  paths  through  a  plausible  bell-shaped  nozzle.  The  6 

upper  two  show  rays  that  penetrate  the  throat;  the  lower  two  show  rays  that  do 

not.  (The  last  one,  lower  right,  shows  a  "‘whispering  gallery"  ray  path,  highly 
atypical.) 

4  Gain  factors  for  bell-shaped  and  conical  nozzles.  (The  annotations  are  the  6 

coefficients  of  the  polynomials  defining  the  nozzle  shapes.) 

5  The  nozzle  shapes,  together  with  their  defining  polynomial  coefficients,  used  for 
generating  the  ray-trace  images  of  Figure  3  and  the  nozzle  gain-factor  curves  of 
Figure  4. 

6  The  on-axis  gain  factor  for  a  "long"  conical  nozzle.  The  total  on-axis  RCS  8 

enhancement  due  to  the  nozzle  is  twice  this.  For  a  conical  nozzle  of  arbitrary 

length,  the  on-axis  gain  factor  is  the  lesser  of  the  long-nozzle  gain  factor  and  the 
mouth-to-throat  area  ratio. 

7  The  extended-return  decay  rates  of  various  liquid-fueled  rocket  engines,  9 

domestic  and  foreign,  inferred  from  their  characteristic  chamber  lengths. 

8  The  energy  captured  by  the  combustion  chamber  is  the  noncoherent  integral  13 

over  the  hemispherical  surface  C  of  the  energy  flux  crossing  it. 


IX 


9 


15 


Throat  penetration  of  a  conical  nozzle  for  axial  incidence.  All  rays  within  the 
radius  6 'penetrate  the  nozzle  throat.  The  direct  rays  (blue)  lie  within  the  central 
circle  of  radius  a\  the  single-reflection  rays  (red)  lie  outside  that  circle  in  the 
annulus  of  upper  radius  b'  If  the  cone  half  angle  is  more  than  22.5°,  only  those 
two  ray  types  can  penetrate  because  the  second  reflection  sends  any  ray 
obliquely  back  toward  the  nozzle  mouth. 

10  If  the  nozzle  cone  half  angle  is  less  than  22.5°,  double-reflection  rays  can  also  16 

penetrate  the  throat.  These  rays  (green)  lie  in  a  second  annulus  lying  outside  the 
one  occupied  by  the  single-reflection  rays  (red).  After  the  second  reflection, 
they  are  still  heading,  obliquely,  toward  the  nozzle  throat  Triple-reflection  rays 
achieve  throat  penetration  only  for  cone  half  angles  smaller  than  15°. 


x 


1.  CAVITY  SIGNATURE:  A  THERMAL  MODEL 


1.1  INTRODUCTION 

Cavities  that  generate  long  delayed  radar  returns  present  a  difficult  problem  for  radar  signature 
modeling.  Physically  such  cavities  are  either  complicated  (the  end  of  a  separated  attitude  control  module 
with  a  crowded  interior,  for  example)  or  they  have  the  combination  of  a  small  aperture  and  a  large  volume 
(the  combustion  chamber  of  a  liquid-fuel  rocket  or  the  burned-out  case  of  a  solid-fuel  rocket).  The  long 
delay  derives,  in  the  first  case,  from  the  multiple  scattering  from  objects  within  the  cavity,  and  in  the 
second,  from  the  slow  leakage  imposed  by  the  small  aperture.  Modeling  such  cavities 
electromagnetically,  therefore,  is  simultaneously  more  important  than  modeling  simple  cavities  but  also 
more  difficult 

A  recent  review  by  Anastassiu1  of  the  various  direct  approaches  to  modeling  cavities  concludes  that 
they  all  fail  when  the  complexity  of  the  cavity  reaches  a  certain  level.  He  describes  exact  methods  using 
differential  equations  and  integral  equations,  modal  techniques,  finite  element  methods,  and  the  high- 
frequency,  spectral,  and  ray-tracing  methods.  Unfortunately,  the  level  of  complexity  at  which  they  fail  is 
still  short  of  the  levels  reached  by  cavities  of  practical  interest. 

The  thermodynamic  model  turns  away  from  all  these  direct  approaches.  It  approaches  the  problem 
from  the  other  direction,  by  taking  the  complexity  and  slow  leakage  to  their  limits.  Specifically,  it 
introduces  the  concept  of  a  well-mixed  thermal  cavity  within  which  the  captured  radar  energy  is 
distributed  uniformly  in  space,  in  direction  of  flow,  and  in  polarization.  It  also  assumes  that  the  cavity  is 
lossless.  And,  finally,  it  assumes  that  all  dimensions  of  the  cavity  are  large  compared  to  a  wavelength. 

The  consequences  of  these  assumptions  are  simple  formulas  for  the  cavity’s  radar  signature, 
narrowband  and  wideband.  They  apply  to  the  basic  well-mixed  cavity  for  which  the  plane  of  the  cavity 
aperture  has  full  hemispheric  exposure.  The  formulas  depend  only  on  the  aperture  area,  the  cavity 
volume,  and  the  aspect  angles,  incident  and  scattering. 

A  rocket  nozzle  modifies  the  simple  model.  It  shields  the  aperture  (the  nozzle  throat)  from  the 
incoming  radar  energy  at  large  aspect  angles  and  increases  its  effective  energy-collecting  area  at  small 
aspect  angles.  It  has  the  reciprocal  effect  on  the  energy  discharge  from  the  cavity,  in  that  it  narrows  and 
intensifies  its  distribution  In  fact,  as  will  be  demonstrated,  the  same  numerical  function  of  aspect,  called 
here  the  nozzle  gain  factor,  modifies  both  the  energy  collection  and  energy  discharge. 

An  example,  applying  the  model  to  a  fictitious  rocket  engine  with  a  bell-shaped  nozzle,  shows  how 
large  the  magnifying  effect  of  the  nozzle  can  be.  With  the  unexceptional  dimensions  of  a  20°  cone  half 


1  H.  T.  Anastassiu,  “A  Review  of  Eleetromagnetie  Scattering  Analysis  for  Inlets,  Cavities,  and  Open  Duets,"  IEEE 
Antennas  and  Propagation  Magazine,  Vol.  45,  No.  6,  pp.  27-40  (Dee.  2003). 


1 


angle  and  a  9: 1  area  ratio  of  the  mouth  to  the  throat,  the  nozzle  increases  the  on-axis  backscattering  RCS 
of  the  combustion  chamber  by  some  1 8  dB. 

1.2  BASIC  THERMAL  CAVITY 

Figure  1  shows  the  generic  basic  cavity  of  volume  V  with  a  plane  aperture  of  area  A.  If  the  incident 
energy  flux  density  (Joules/m2)  is  Wj ,  and  the  aperture  dimensions  are  large  compared  with  a  wavelength, 
then  the  total  energy  passing  through  the  aperture  into  the  cavity  is  simply  fVA cos#  .  ( 6 t  is  the  angle 
between  the  transmitting  radar’s  line  of  sight  and  the  aperture-plane  normal.)  If  the  cavity  is  lossless,  all 
of  this  energy  leaks  back  out  again,  its  rate  of  leakage  determined  by  the  internal  energy  density  and  the 
aperture  area,  taking  account  of  all  the  different  directions  with  which  it  approaches  the  aperture.  Section 
A.l  does  the  algebra,  leading  to  the  follow  ing  simple  formula,  independent  of  polarization  and  frequency, 
for  the  narrowband  bistatic  RCS  of  the  cavity: 


rr  =  2/tcos#  cos#  , 

where  6S  is  the  angle  between  the  receiving  radar’s  line  of  sight  and  the  aperture-plane  normal. 


Figure  l.  The  basic  cavity ,  showing  the  aperture,  the  single  incidence  direction,  and  the  energy  discharging  along 
multiple  scattering  directions. 


The  frequency  independence  is  based  on  the  assumption  that  all  dimensions  of  the  cavity  and  the 
aperture  are  large  compared  with  a  wavelength.  The  smallest  dimension,  and  therefore  the  limiting 
dimension,  is  likely  to  be  that  of  the  aperture.  If  it  is  not  large,  in  wavelengths,  then  its  effective  area 
would  be  different  from  its  actual  area.  And  if  it  is  below  cut-off  for  electromagnetic  wave  transmission, 
no  energy  would  enter  the  combustion  chamber,  and  there  would  be  no  extended  return  at  all.  Fortunately, 
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since  the  cut-off  radius  is  given2 *  very  closely  by  ka  =  1 .8  (k  is  the  wave  number  and  a  the  throat  radius), 
at  10  GHz  (X-band)  the  eut-off  diameter  2 a  is  1.7  em.  This  is  small  compared  to  the  throat  diameter  of 
large  roeket  motors,  so  for  the  higher-frequency  radars  we  ean  expect  the  “large  aperture"  assumption 
usually  to  be  satisfied.  Moreover,  the  eut-off  is  very  abrupt \  Above  eut-off,  the  transmission  is  essentially 
perfeet;  below  eut-off,  there  is  none. 

The  wideband  radar  return  from  the  cavity  exhibits  an  extended  return  stretehing  out  down-range  of 
the  eavity  aperture.  The  measured  length  of  the  return  is  controlled  by  its  deeay  rate,  and  that,  in  turn,  is 
determined  by  the  rate  of  leakage  out  through  the  aperture  of  the  energy  initially  eaptured  by  the  cavity. 
Section  A. 2,  starting  from  the  assumption  that  the  incident  radar  energy  instantaneously  fills  the  cavity 
with  a  well-mixed  thermal  radiation  field,  derives  this  rate  of  leakage  and  then,  by  balancing  the  input  and 
output  energies,  obtains  the  following  explicit  formula,  plotted  in  Figure  2,  for  the  wideband  signature  of 
the  cavity's  extended  return: 


Act  = 


A rA2  cos 0t  cos 0S 
V 


c\\)(-Ar!2V)  - 


Here,  Fis  the  cavity  volume,  A  is  the  aperture  area,  0.  and  6s  are,  respectively,  the  angles  between 
the  aperture-plane  normal  and  the  radar  lines  of  sight,  transmitting  and  receiving,  Ar  is  the  radar's  range 
resolution,  and  r  is  the  down-range  distance  from  the  cavity  aperture.  The  ratio  V/A.  which  occurs  twice  in 
this  formula,  is  the  quantity  defined  by  Sutton  and  Biblarz4  in  their  treatise  on  rocket  propulsion  as  the 
characteristic  chamber  length  L  of  the  cav  ity. 

In  Figure  2,  the  line  plot  (labeled  “original  formula")  of  this  expression  is  normalized  by  A r, 
applies  to  on-axis  baekseattering,  and  uses  the  unexceptional  values  0.2  m  for  the  equivalent  radius  of  the 
cavity  and  0.1  m  for  the  aperture  radius.  Its  unrealistic  abrupt  start  is  the  result  of  the  assumption  that  the 
well-mixed  energy  in  the  eavity  is  established  instantaneously.  The  figure  also  shows  a  second  line  plot 
that  includes  a  proposed  eorreetion  for  this  problem. 

The  correction  takes  account  of  the  fact  that,  initially,  the  bulk  of  the  energy  in  the  cavity  is  mov  ing 
internally  away  from  the  aperture,  and  so  little  is  flowing  baek  out.  This  has  the  effect  of  delaying  the 
apparent  range  at  whieh  the  peak  of  the  eav  ity  return  occurs.  The  rationale  for  the  form  of  the  adjustment 
is  that  if  the  eav  ity  is  the  interior  of  a  simple  wide-open  tin  can,  the  range  delay,  for  normal  incidence,  is 
simply  the  depth  of  the  can.  This  depth  is  also  equal  to  the  volume  of  the  can  V  divided  by  the  area  A  of 
the  aperture.  That  is,  the  depth  is  the  characteristic  chamber  length  L* .  This  suggests  that  the  formulas 
above  for  Act  should  be  multiplied  by  the  factor  1  -  cxp(-r  /  C ) ,  whieh  rises  from  zero  at  r  =  0  and 
approaches  one  for  large  r.  This  ehanges  the  formula  for  the  distributed  RCS  of  the  cav  ity  to 


2  A.  Roberts,  “Electromagnetic  Theory  of  Diffraction  by  a  Circular  Aperture  in  a  Thick,  Perfectly  Conducting 
Screen,”  J.  Opt .  Soc  Am..  Vol.  4,  No  10,  Oct.  1987. 

1  Ibid. 

4  G.  P.  Sutton  and  O.  Biblarz,  Rocket  Propulsion  Elements ,  John  Wiley  &  Sons,  New  York,  2001 . 
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Figure  2.  The  wideband  art-axis  backse ottered  cavity  signature  for  the  basic  well-mixed  thermal  cavity  The 
unrealistic  abrupt  start  of  the  original  formula  is  suppressed  in  the  modified  formula.  The  \ -axis  units  are  square 
meters  of  radar  cross  section  per  meter  of  range,  expressed  in  dB  (The  effective  cavity  radius  is  the  radius  of  a 
sphere  having  the  same  volume  as  the  cavity.) 


r 

Act  =  :^pr  A  cos 0  cos#  [1  -  exp (-r !  C )]  exp(-r  /2L* ) , 

in  which  now  every  occurrence  of  the  ratio  V/A  has  been  replaced  by  l' .  The  factor  3/2  ensures  that  the 
narrowband  RCS  remains  at  cr  =  2  A  cos #  cos#  . 

The  asymptotic  decay  rate  of  the  extended  return  in  dB/m  is  given  by  the  coefficient  of  v  in  the 
expression  lOlog10A(T;  it  is  (5/ L*)logl0e  or,  equivalently,  2  17/  L*  dB/m  For  the  cavity  and  aperture 
dimensions  used  in  Figure  2,  the  rate  is  2.03  dB/m. 

Another  potentially  significant  adjustment  to  the  formulas  accommodates  the  more  general 
situation  in  which  the  cavity  is  not  lossless.  Specifically,  the  additional  exponential  factor  exp(-f  /  r)  is 
appended  to  the  wideband  formulas  above,  where  T  is  the  energy-decay  time  constant  the  cavity  would 
have  if  closed.  That  is,  the  exponential  factor  e\p(-r  /  2L* )  is  replaced  by  cxp(-r  /  2ll  -  2r  /  ct)  ,  w  here 
c  is  the  speed  of  light,  and  the  equivalence  relationship  between  time  and  range  has  also  been  applied.  The 
corresponding  adjustment  to  the  narrowband  formula  is  to  multiply  the  lossless  RCS  by  the  factor 
(1  +  4  L*  fer)  1  . 
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1.3  NOZZLE  GAIN  FACTOR 


A  gain  factor  Ft{6t)  applied  to  the  formula  for  the  RCS  of  the  basic  cavity  defines  the  way  the 
nozzle  of  the  rocket  engine  reshapes  the  energy-gathering  pattern  of  the  basic  cavity  aperture.  A  second 
gain  factor  Fs(0s)  defines  the  reshaping  of  its  energy  discharge.  Section  A. 3  shows  that,  by  reciprocity, 
these  factors  are  identical  the  subscripts  /'  and  5  can  be  dropped.  The  narrowband  and  wideband  bistatic 
RCS  formulas  for  the  basic  thermal  cavity  thus  become 


a  =  2 AF{0  )F(0s) cos#  cos6> 

and 

Act  =  AF(dl  )F(0s) cos 6?  cos&s  (l  -  exp (-r/L*  ))cxp(-r/2/Z )  , 
respectively,  when  applied  to  the  rocket  engine. 

Ray  tracing  can  be  used  to  evaluate  the  gain  factor  by  counting  how  many  more  rays,  of  a  uniform 
parallel  flux  of  incident  rays,  pass  through  the  aperture  with  the  nozzle  present.  It  involves  tracing  the  ray 
reflections  off  the  nozzle  walls.  Section  A. 5  illustrates  the  process  for  a  conical  nozzle.  Section  A. 6  is  the 
MATLAB  function  for  evaluating  the  gain  factor  of  a  nozzle  shape  defined  by  an  arbitrary  polynomial. 

Figure  3  shows  the  result  of  applying  this  MATLAB  function  to  tracing  four  different  rays 
interacting  with  a  plausible  bell-shaped  nozzle.  The  rays  in  the  upper  two  images  penetrate  the  throat 
directly  or  after  two  reflections;  those  in  the  lower  two  reflect  back  out  without  penetrating  the  throat. 
(The  trajectory  of  the  highly  atypical  ray  in  the  lower  right  image  propels  it  along  a  “whispering  gallery” 
path,  undergoing  multiple  reflections  as  it  moves  around  the  inside  wall.) 

Two  gain  factor  curves  computed  by  the  function,  one  for  the  bell-shaped  nozzle  and  one  for  a 
similar  conical  nozzle,  are  presented  in  Figure  4  They  show  that  a  nozzle  can  have  a  substantial  one-way 
on-axis  gain  in  this  case  nearly  10  dB.  The  total  nozzle  gain  is  twice  this.  The  nozzle  achieves  this  gain 
by  narrowing  the  angular  field  both  from  which  energy  is  gathered  and  to  which  it  is  discharged.  These 
two  gain-factor  curves  make  the  important  point  that  the  difference  between  the  gain  factors  of  a  bell- 
shaped  nozzle  and  of  its  conical  nozzle  approximation  can  be  small  enough  that  the  latter  will  often  be 
accurate  enough. 

Figure  5  shows  the  nozzle  shapes  themselves.  The  area  ratio,  mouth  to  throat,  of  both  nozzles  is  9, 
and  the  cone  half  angle  of  the  conical  nozzle  is  20°. 
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Figure  3.  Ray  paths  through  a  plausible  bell-shaped  nozzle.  The  upper  two  show  rays  that  penetrate  the  throat;  the 
lower  two  show  rays  that  do  not .  (The  last  one,  lower  right  shows  a  "whispering  gallery  "  ray  path,  highly  atypical.) 


Figure  4.  Gain  factors  for  bell-shaped  and  conical  nozzles.  (The  annotations  are  the  coefficients,  m  MATLAB 
format,  of  the  polynomials  defining  the  nozzle  shapes.) 
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Figure  5.  The  nozzle  shapes,  together  with  their  defining  polynomial  coefficients,  used  for  generating  the  ray-trace 
images  of  Figure  3  and  the  nozzle  gain-factor  curves  of  Figure  4. 


A  useful  measure  of  the  gain  factor  curve  is  its  angular  half  width  0W  .  It  indicates  the  extent  to 
which  the  nozzle  narrows  the  view  of  the  cavity,  and  it  is  readily  estimated  from  the  on-axis  value  of  the 
gain  factor.  Section  A. 4,  relying  on  the  principle  of  energy  conservation,  demonstrates  that  the 
relationship  is: 


9W  =arcsin(F  1  2(0)). 

But  the  gain  factor  of  the  actual  nozzle  can  he  evaluated  with  reasonable  accuracy  by  considering 
its  conical  nozzle  approximation.  In  particular,  as  the  next  section  describes,  the  on-axis  gain  factor  of  a 
conical  nozzle  can  be  expressed  in  simple  closed  form.  That  provides  a  method  for  evaluating  0W  without 
any  numerical  ray  tracing. 

For  the  case  represented  by  the  bell-shaped  nozzle  of  Figure  4,  this  half  width  works  out  to  be 
19.5°,  which  corresponds  to  a  point  on  the  gain-factor  curve  4.2  dB  below  its  peak  value.  For  the 
backscattering  RCS,  which  involves  the  square  of  the  gain  factor  as  well  as  the  square  of  cos  0  ,  this  point 
is  9.0  dB  below  the  peak. 

1.4  ON-AXIS  GAIN  FACTOR  FOR  A  CONICAL  NOZZLE 

The  situation  is  particularly  simple  for  a  “short”  purely  conical  nozzle  with  a  cone  half  angle  of  less 
than  45°.  At  axial  incidence,  every  ray  incident  within  the  nozzle  mouth  penetrates  the  throat  and 
soF(0)  =  b2/a2 ,  where  a  and  b  are,  respectively,  the  throat  and  mouth  radii.  But  if  the  nozzle  is  “long,” 
rays  outside  a  certain  radius  get  reflected  back  out  of  the  nozzle  and  so  do  not  penetrate  the  throat.  The 
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effective  energy-collection  radius  b'  of' the  nozzle  is  then  less  than  h,  and  F(0)  =  b'~/a  In  general, 
therefore. 


F( 0)  =  mii r(b,b')fa2 . 


Section  A. 5  shows  that  for  the  "long”  nozzle  with  a  half-angle  y  between  22.5°  and  45°, 
A//a  =  4cos2  y-l ,  and  if  y  lies  between  15°  and  22.5°,  b'ia  =  Ibeos4  y—  l2cos:  y+  1  .  Extending  this 
process  to  even  smaller  eone  half  angles  produees  the  plot  in  Figure  6. 


Figure  6.  The  on-axis  gain  factor  for  a  "long”  conical  nozzle.  The  total  on-axis  RCS  enhancement  due  to  the  nozzle 
is  twice  this .  For  a  conical  nozzle  of  arbitrary  length ,  the  on-axis  gain  factor  is  the  lesser  of  the  long-nozzle  gain 
factor  and  the  mouth-to-throat  area  ratio. 


The  graph  shows  that  the  nozzle  can  magnify  the  axial  backscattering  RCS  of  the  combustion 
chamber  by  a  large  factor — some  18  dB  in  the  case  of  a  20°  cone  half  angle,  if  the  radius  ratio,  mouth  to 
throat,  of  the  nozzle  is  no  smaller  than  7.94.  (The  RCS  enhancement  is,  in  dB,  twice  the  gain  faetor.) 

1.5  DISCUSSION 

The  ideal  thermal  model  described  above  is  limited.  It  cannot  model  simple  cavities.  For  example, 
the  narrowband  RCS  of  a  simple  can-shaped  cavity  at  normal  incidence  is  about  4  ttA2/A2  (the  RCS  of  its 
base),  whereas  the  corresponding  thermal  RCS  is  2 A.  At  X-band,  if  the  diameter  is  1  m,  these  formulas 
imply  the  very  different  RCS  values  of  45.4  and  2.0  dBsm,  respectively.  The  high  number  is  not 
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significantly  reduced  if  we  include  a  small  internal  scatterer — much  more  complexity  is  necessary  before 
the  ideal  thermal  model  becomes  adequate. 

It  should  also  be  noted  that  the  cavity  need  not  be  complex  to  generate  an  extended  return. 
Dumaman  et  al.  show  a  simple  cylindrical  cavity  with  an  extended  return  stretching  out  behind  the  target 
more  than  fifteen  times  its  physical  depth,  if  the  aspect  angle  is  large.  And  they  show  that  the  extended 
return  is  not  present  at  smaller  aspect  angles-  very  nonthermal  behavior.  Even  though  the  energy  within 
the  cavity  undergoes  multiple  bounces,  the  simple  geometry  of  the  cavity  ensures  that  the  bounces  are 
very  ordered.  The  energy  is  not  well  mixed. 

On  the  other  hand,  if  it  is  true  that  the  energy  is  well  mixed,  then  all  that  follows  in  producing  the 
RCS  formulas  would  seem  to  be  logically  necessary.  Clearly  wanting  here  are  comparisons  of  the  theory 
with  some  measured  cavity-return  data,  but  that  must  for  the  moment  be  handled  in  another  context. 

However,  some  testable  predictions  of  the  decay  rate  can  be  made  from  the  values  of  L  published 
for  various  real  liquid-fueled  rocket  engines,  domestic  and  foreign,  by  Sutton  and  Biblarz* * 6.  They  arc 
displayed  in  the  graph  in  Figure  7. 


Figure  7.  The  extended-return  decay  rates  of  various  liquid-fueled  rocket  engines,  domestic  and  foreign .  inferred 
from  their  characteristic  chamber  lengths . 


A.  J  Dumanian,  E.  C.  Burt,  and  B.  A.  Kemp,  “A  Component  Model  Approach  for  the  RCS  Validation  of  an 

Electrically  Large  Open-Ended  Cavity"  (to  be  published). 

6  G.  P.  Sutton  and  C.  Biblarz,  op.  c/7.,  p.  272,  392. 
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A  remarkable  feature  of  these  numbers  is  the  relatively  narrow  range  over  whieh  the  extended- 
return  decay  rates  of  the  high-thrust  engines  lie.  These  six  engines,  the  largest  of  whieh  has  a  thrust  more 
than  180  times  that  of  the  smallest,  differ  in  their  predicted  decay  rates  by  a  ratio  no  more  than  1 .64.  How 
well  these  theoretical  decay  rates  agree  with  the  aetual  ones  is  yet  to  be  determined. 

A  power-law  curve  with  the  engine's  vaeuum  thrust  as  the  independent  variable  is  shown  fitted  to 
the  plotted  points.  If  nothing  else,  it  provides  a  framework  for  the  numbers. 

1.6  CONCLUSIONS 

The  well-mixed  eavity  model  produces  simple  formulas  for  the  extended  radar  signatures  of  liquid- 
fueled  roeket  engines  and  highly  cluttered  cavities.  The  energy  capture  and  energy  release  processes  of 
the  cavity  are  so  simple  to  model  that  the  overall  model  accuracy  depends  wholly  on  the  assumption  that 
the  eavity  energy  is  well  mixed  a  seemingly  safe  assumption  for  these  cavities.  However,  this  report 
includes  no  comparison  with  measured  data. 
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APPENDIX 


A.l  BASIC  WELL-MIXED  THERMAL  CAVITY 

The  energy  E  captured  by  the  basic  cavity  is  given  by  E  =  IV.  A  cos  0t ,  where  Wt  is  the  incident 
energy  flux  density,  A  is  the  area  of  the  cavity's  aperture  and  0i  is  the  incidence  angle  (the  angle  between 
the  line  of  sight  to  the  transmitting  radar  and  the  normal  to  the  plane  of  the  aperture). 

The  perfect  mixing  assumption  for  the  energy  within  the  cavity  implies  that  the  energy  within  is 
uniformly  distributed  in  its  direction  of  propagation.  That  in  turn  implies  that  the  reradiated  energy  flux 
density  \Vs  in  any  direction  is  proportional  to  the  projected  area  of  the  aperture  in  that  direction.  In  other 
words,  it  is  reradiated  with  a  cos#v  dependence  on  the  scattering  angle  6?  (the  angle  between  the  line  of 
sight  to  the  receiving  radar  and  the  normal  to  the  plane  of  the  aperture),  namely  Ws  =a cos#  ,  with  a  to 
be  determined. 

The  total  energy  reradiated,  for  a  lossless  cavity,  can  be  equated  to  the  energy  captured.  That  is,  the 

reradiated  energy  flux  density  Ws  must  satisfy  the  equation  E  =  R2  j \\Vs(Kls  .  Carrying  out  this  integral 
over  the  hemisphere,  with  Ws  represented  as  tfcos^J. ,  allows  a  to  be  evaluated,  and  IV  to  be  expressed 
as  Wx  =  fVjAcosfy  cos#  /  kR  .  This,  together  with  the  definition  Ws  -  Wp  l(AnRr)  for  the  radar  cross 

section,  yields  the  required  expression  for  the  bistatic  narrowband  total  radar  cross  section  of  the  basic 
well-mixed  thermal  cavity:  <7  =  44  cos#  cos#  . 

Important  to  note  here  is  the  fact  that  this  expression  represents  the  total  energy  scattered  by  the 
cavity.  According  to  the  perfect  mixing  assumption,  this  must  be  divided  equally,  statistically,  between 
whichever  two  orthogonal  polarizations  one  chooses  to  define.  The  formula  representing  the  true 
measured  narrowband  cross  section  therefore,  evaluated  at  any  polarization,  is  just  half  this,  namely 

<7  =  2  A  cos  0i  cos#. . 


A. 2  WIDEBAND  SIGNATURE 

Since  the  cavity  is  assumed  to  be  lossless,  the  rate  of  decay  of  the  energy  in  the  cavity  is  equal  to 
the  total  power  flux  back  out  through  the  aperture.  This  can  be  evaluated  from  the  total  energy  E  currently 
in  the  cavity  using  the  well-mixed  cavity  assumption. 

Specifically,  the  directional  energy  density  in  the  cavity  is  E  /  4rtV  Joules/m3/steradian,  where  Lis 
the  cavity  volume,  so  the  directional  power  flux  density  is  Ec  !  AkV  watts/m2/steradian,  where  c  is  the 
speed  of  light,  and  the  directional  power  flux  thorough  the  aperture  is  EcA cos 0  /  4kV  watts/steradian, 
where  A  is  the  aperture  area  and  0  is  the  angle  between  the  propagation  direction  and  the  aperture-plane 
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normal.  The  total  power  flux  through  the  aperture  is  the  integral  of  this  over  the  external  hemisphere,  that 
is,  (cA  /  4V)E  watts.  This  is  the  negative  of  the  rate  of  change  of  the  total  energy  in  the  cavity,  implying 
that  the  total  energy  in  the  cavity  decays  according  to  the  law 

E  —  W.A cos#  exp(-c/fC  /  4V) , 

where  the  factor  combination  W.A cos#  represents  the  initial  energy  captured  by  the  cavity. 

The  result  of  substituting  this  expression  for  E  into  the  formula  EcAcosO  / 4nV  from  the  previous 
paragraph  for  the  directional  power  flux  thorough  the  aperture  in  the  direction  of  the  receiving  radar 
converts  it  into 


cvCffcos#,  cos#v 
4nV 


exp(-cyff ! 4V)  • 


Its  units  are  watts/steradian. 

Multiplying  this  same  expression  by  the  small  time  increment  At  (the  range  resolution 
Ar  expressed  in  units  of  time)  gives  the  combined  incremental  radar-directed  energy  flux 
(Joules/steradian)  of  both  polarizations,  which,  from  the  definition  of  radar  cross  section,  is  also  equal  to 
twice  the  incremental  radar  cross  section  Act  measured  in  any  one  polarization  times  the  incident  energy 
flux  W  (Joules/m  )  divided  by  Air.  Solving  this  equation  for  Act,  and  observing  the  equivalence  of  ct 
with  2 1\  finally  yields  the  following  expression  for  the  bistatic  incremental  wideband  RCS  of  the  basic 
well-mixed  cavity: 


ArA  cos#  cos# 

A  a  = - - -cxp(-Ar/2V)  • 


The  narrowband  RCS  is  the  range  integral  of  Act /A  /  ,  namely  cr  =  2,4  cos  #r  cos#  ,  which  agrees,  as 
it  should,  with  the  direct  derivation  of  the  previous  section. 

Here,  V  is  the  cavity  volume,  A  is  the  aperture  area,  0t  and  #  are,  respectively,  the  angles  between 
the  aperture-plane  normal  and  the  radar  lines  of  sight,  transmitting  and  receiving,  and  A)'  is  the  radar’s 
range  resolution. 

The  asymptotic  decay  rate  of  the  extended  return  in  dB  m  is  given  by  the  coefficient  of  r  in  the 
expression  101og)n  Act  ;  it  is  (5/1  /  V)\og]0e  or,  equivalently,  2. 1 1A !  V  dB/m. 
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A.3  NOZZLE  GAIN-FACTOR  RECIPROCITY 


The  reciprocity  theorem  implies  that  the  nozzle  gain  factor  shaping  a  combustion  chamber’s  energy 
capture  is  the  same  function  of  aspect  as  the  nozzle  gain  factor  shaping  the  combustion  chamber's  energy 
discharge. 

For  proof  of  this  statement,  it  is  useful  to  imagine  the  nozzle  alone,  unattached  to  its  combustion 
chamber.  Then,  the  radar-pulse  energy  captured  by  the  cavity  is  the  noncoherent  integral  over  the 
hemispherical  surface  C  of  the  energy  penetrating  the  nozzle  throat  (see  Figure  8).  In  the  other  direction, 
the  total  energy  radiating  from  the  mouth  of  the  cavity  back  in  the  direction  of  the  radar  is  the 
noncoherent  integral  over  the  hemispherical  surface  C  of  the  energy  radiated  in  the  direction  of  the  radar 
by  each  one  of  a  uniform  distribution  of  identical  emitters.  (The  emitters  are  identical  and  uniformly 
distributed  over  C  because,  according  to  the  assumptions  of  the  thermal  cavity  model,  the  radiation  field 
inside  the  cavity  is  uniformly  distributed  in  direction.) 


Figure  8.  The  energy’  captured  by  the  combustion  chamber  is  the  noncoherent  integral  over  the  hemispherical 
surface  C  of  the  energy  flux  crossing  it. 


The  reciprocity  theorem  equates,  on  the  one  hand,  the  field  at  each  point  on  surface  C  due  to  the 
radar,  and,  on  the  other,  the  field  at  the  radar  due  to  each  point  emitter  on  surface  C.  This  equality 
extends,  therefore,  to  the  two  integrals. 

A.4  RCS  LOBE-WIDTH  ESTIMATE 

The  expression  F(6?)cos#  involving  the  nozzle  gain  factor  fully  describes  the  directional 
properties  of  the  engine’s  radar  signature.  But  a  rough  indication  of  the  main-lobe  half  width  0W  of  the 
RCS  can  be  derived  that  depends  only  on  the  on-axis  value  F(0) .  The  derivation  depends  on  the 

statement  E  =  2R2  \ws  d£ls  of  the  conservation  of  energy,  equating  the  total  energy  captured  by  the 
combustion  chamber  to  the  total  energy  reradiated.  Here  2 Ws  is  the  total  reradiated  energy  flux  density 
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( IV  in  each  polarization)  as  a  function  of  solid  angle  ils  and  the  integration  is  carried  out  over  the 
hemispherical  surface  of  (large)  radius  R  across  which  the  energy  passes. 

The  application  to  this  conservation  statement  of  the  three  identities  E  =  IV  AF(9) cos9  (for  the 
total  energy  captured  by  the  combustion  chamber),  Ws  a  l(4/r R  )  (for  the  definition  of  radar  cross 
section),  and  a  =  2AF(9,)F(9)cos9tcos9  (for  the  radar  cross  section  of  the  combined  combustion 
chamber  and  nozzle  for  each  polarization),  reduces  it  to  1  =/r  jV(#  )cos#c/£2v .  (A  is  the  nozzle  throat 
area  and  IV  the  incident  energy  flux  density.)  And  since  dQs  =  sin  9scWsd(ps ,  this  can  be  rewritten 

It  2  &h 

I  =  jV(0,)sin  2 etdes  =  F(0)  Jsin  2 0td0s  - 

0  0 

where  9W  is  the  main-lobe  half  width.  Evaluating  the  integral  and  then  solving  for  9W  leads  to  the  result 

9h  =  arcsin(F  1 2  (())). 

A.5  ON-AXIS  GAIN  FACTOR  FOR  A  CONICAL  NOZZLE 

For  a  conical  model  of  the  nozzle,  the  simple  ray  geometry  for  axial  incidence  allows  easy 
evaluation  of  the  gain  factor  on-axis.  If  the  cone  half  angle  y  is  greater  than  45°,  the  nozzle  gain  factor  at 
axial  incidence  is  1 ;  in  other  words,  at  axial  incidence  no  ray  can  penetrate  the  nozzle  throat  after 
reflection  off  the  nozzle  wall  and  so  the  nozzle  has  no  effect.  The  only  energy  to  penetrate  the  nozzle 
throat  does  so  directly,  just  as  it  would  were  no  nozzle  present. 

For  cone  half  angles  less  than  45°,  but  no  smaller  than  22.5°,  the  only  rays  that  penetrate  the  throat 
at  axial  incidence  are  those  that  pass  either  directly  through  the  throat  or  through  it  after  one  reflection  off 
the  nozzle  wall.  The  circular  symmetry  of  the  nozzle  ensures  that  this  additional  collection  area  is  a  circle 
with  a  radius  b'  that  is  bigger  than  the  radius  a  of  the  nozzle  throat.  Every  incident  ray  that  is  within  a 
radius  b'  of  the  nozzle  axis  penetrates  the  throat,  either  directly  or  after  one  reflection  off  the  nozzle  wall. 
Thus,  the  gain  factor  is  simply  the  ratio  of  the  two  circle  areas,  specifically  (b'  /  a)2 .  This  ray  geometry  is 
shown  in  Figure  9.  However,  if  the  nozzle  is  short  enough,  the  radius  of  its  energy  collection  circle  will  be 
just/),  the  radius  of  the  nozzle  mouth.  Thus,  in  general,  the  gain  factor  will  be  {min ty\b)!a)2. 
Evaluating  b'  is  a  simple  exercise  in  trigonometry.  The  final  solution  for  the  gain  factor  for  axial 
incidence  is  then 


F(0)  =  miri[4cos2(y)-lb/a]  ,  (22.5°  <  y  <  45°) , 
where  b  and  a  are  the  mouth  and  throat  radii  of  the  nozzle  and  y  is  the  cone  half  angle. 
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Figure  9.  Throat  penetration  of  a  conical  nozzle  for  axial  incidence.  All  rays  within  the  radius  b'  penetrate  the 
nozzle  throat.  The  direct  rays  (blue)  lie  within  the  central  circle  of  radius  a;  the  single-reflection  rays  (red)  lie 
outside  that  circle  in  the  annulus  of  upper  radius  b '  If  the  cone  half  angle  is  more  than  22.5°,  only  those  two  ray 
types  can  penetrate  because  the  second  reflection  sends  any  ray  obliquely  back  toward  the  nozzle  mouth. 


For  smaller  cone  half  angles,  down  to  15°,  at  axial  incidence,  some  rays  also  penetrate  the  throat 
after  two  reflections  off  the  nozzle  wall,  as  shown  in  Figure  10.  In  this  situation,  following  the  same  logic 
as  for  the  onc-reflection  case,  the  gain  factor  for  axial  incidence  is  given  by 

F(0)  =  mnr  [1 6  cos4  (/)  - 1 2  cos2  ( y)  + 1, b  l a]  ( 1 5°  <  y  <  22.5°) . 

Since  the  total  gain  factor  is  the  product  of  the  collecting  and  discharge  gain  factors,  and  since  these 
factors  are  equal,  the  total  gain  factor  of  the  nozzle  for  backscattering  at  axial  incidence  is  simply  F “(0) . 

This  analysis  is  readily  extended  to  smaller  cone  half  angles  with  their  penetrating  rays  of  even 
greater  bounce  count.  The  resulting  gain  factor  is  presented  graphically  in  Figure  6  in  the  main  text,  in 
which  the  “long  nozzle”  assumption  is  made.  That  is,  it  is  assumed  that  the  nozzle  mouth  radius  is  larger 
than  the  geometric  collecting  radius  determined  by  the  formulas  in  cos(/).  The  graph  shows  that  the 
nozzle  can  magnify  the  axial  backscattering  RCS  of  the  combustion  chamber  by  a  large  factor — some 
1 8  dB  in  the  case  of  a  20°  cone  half  angle,  if  the  radius  ratio,  mouth  to  throat,  of  the  nozzle  is  no  smaller 
than  7.94. 
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Figure  JO.  If  the  nozzle  cone  half  angle  is  less  than  22.5  °f  double-reflection  ravs  can  also  penetrate  the  throat. 
These  rays  (green)  he  in  a  second  anna  his  l\  ing  outside  the  one  occupied  by  the  single-reflation  rays  (red).  After 
the  second  reflection ,  they  are  still  heading ,  obliquely,  tow  ard  the  nozzle  throat.  Triple-reflection  rays  at  hieve 
throat  penetration  only  for  cone  half  angles  smaller  than  15° 


A. 6  RAY-TRACE  MATLAB  FUNCTION 

The  following  MATLAB  function  uses  ray  tracing  to  evaluate  the  gain  factor  at  a  single  angle  of 
incidence  6  for  a  nozzle  whose  shape  is  specified  by  a  polynomial  defining  the  radius  p  of  the  inside 
surface  of  the  nozzle  wall  as  a  function  of  the  axial  dimension  z. 

It  traces  the  paths  of  a  large  number  N  of  parallel  rays  uniformly  distributed,  statistically,  over  the 
nozzle  mouth,  scoring  a  one  for  each  ray  that  eventually  passes  through  the  nozzle  throat  and  a  zero  for 
each  ray  reflected  back  out  of  the  nozzle  mouth  without  ever  passing  through  the  throat.  The  total  score 
Nt  divided  by  N  then  represents  the  fraction  of  the  total  energy  incident  on  the  nozzle  mouth  that  passes 
through  the  throat.  Therefore,  the  nozzle  gain  factor,  defined  as  the  factor  by  which  the  energy  passing 
through  the  throat  with  the  nozzle  present  exceeds  the  energy  passing  through  the  throat  with  the  nozzle 
absent,  is  computed  as  (b1  /  a2)Nf  /  N ,  where  a  and  b  are,  respectively,  the  radii  of  the  nozzle  throat  and 
mouth. 

The  four  images  in  Figure  3  show  the  paths  of  four  rays  computed  using  the  ray-tracing  feature  of 
the  MATLAB  function,  two  of  which  penetrate  the  throat  and  two  that  do  not.  Two  gain  factor  curves 
computed  by  the  function,  one  for  the  bell-shaped  nozzle  and  one  for  a  conical-nozzle  approximation  to  it, 
arc  presented  in  Figure  4. 
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function  G-PolyNozzleGain(c,thdeg,N) 


%  G  is  the  gain  factor  (not  in  dB)  of  the  nozzle  at  the  aspect  angle  thdeg. 

%  The  radius  rho  of  the  nozzle  as  a  function  of  the  axial  dimension  z  is  defined 
%  by  the  polynomial  c:  rho  =  c(l)*zAM  4- ...  +c(M)  *z  +  c(M+l). 

%  The  nozzle  throat  is  at  z  ~  0;  the  nozzle  mouth  is  at  z  =  1 . 

%  N  rays  are  used,  uniformly  distributed  (statistically)  over  the  nozzle  mouth 
%  c=[- 1.7663  5.8506  -7.0398  3.4583  -0.13876  0.18199];  %  9-to-l  bell  nozzle 
%  c=[0. 364  0.182];  %  9-to-l  cone  nozzle 

a~polyval(c,0);,b=poly  val(c,  1 ); 

G^0;,n=0; 
while  n<N 

x=b*(2*rand- 1  );,y=b*(2*rand- 1 ); 
if  xA2+yA2<=bA2 
n=n+l; 

G=G+PolyNozzleTrace(c,thdeg,x,y); 

end 

end 

G=G*bA2/N/aA2; 

function  K~PolyNozzleTrace(cdhdeg,x,y) 

%  checks  to  see  whether  a  particular  ray  incident  on  the  nozzle  mouth 
%  penetrates  throat 

%  c  lists  the  polynomial  coefficients  for  the  nozzle  contour:  rho  =  c(l  )*zAM  +  ...  %  +c(M)  *z  +  c(M+l ) 
%  thdeg  is  the  aspect  angle 

%  x,y  are  the  coordinates  of  the  point  at  which  the  incident  ray 
%  crosses  the  plane  of  the  nozzle’s  exit  aperture 
%  K  is  1  if  the  ray  penetrates  the  throat  and  0  if  it  does  not 

M~length(c)- 1  ;,th~thdeg*pi/l  80; 

A~PolyBase(conv(c,c));  %  the  coefficient  generator  matrix 
rl-[x  y  1];  %  the  starting  point  of  the  ray  in  the  plane  of  the  nozzle  mouth 
pl=[-sin(th)  0  -cos(th)];  %  the  starting  direction  of  the  ray 
for  k-0: 100 
r0=rl;,p0=pl; 

Ls~roots([zeros(l,2*M-2)  l-p0(3).A2  2*r0(l :2)*p0(  1 :2)’  r0(l :2)*r0(l :2)']- ... 

(r0(3)  A(0:2*M))*A.*(p0(3).Afliplr(0:2*M))); 
l^find(imag(Ls): — ^);,Lre=Ls(I);  %  discard  the  complex  roots 
l~find(Lre>le-6);,Lpos=Lre(I);  %  discard  all  non  positive  definite  values  of  L 
L_min(Lpos);  %  retain  only  the  first  encounter  with  the  surface 
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if  length(L):_^0  &  p0(3)<0,  %  no  bounce  point;  ray  passes  directly  through  throat 
K=l;,retum,end 

if  length(L)==:0  &  p0(3)>=0  %  no  bounce  point;  ray  is  reflected  back 
K=0;,retum,end 

rl~rO+L*pO;  %  the  bounce  point 

if  r  1  (3)<0  %  bounce  point  beyond  throat;  ray  passes  directly  through 
K=1  ;,retum,end 

if  r  1  (3)>  1  %  bounce  point  before  mouth;  ray  is  reflected  back 
K=-0;,retum,end 

t=polyval(polyder(c),rl(3));  %  the  surface  slope  at  the  bounce  point 
nl=[rl(  1 :2)  -t*sqrt(rl(l  :2)*rl(  1 :2)*)];  %  the  surface  normal  at  the  bounce  point 
nl_nl/sqrt(nl  *n  1 ');  %  the  unit  surface  normal 
pl~p0-2*nl  *(nl  *p()');  %  the  pointing  direction  of  the  reflected  ray 
end 

%  the  ray  has  exceeded  100  bounces,  so  it's  counted  as  not  passing  through  the  throat 
K-0; 

function  A^PolyBase(PVl) 

%  A  is  the  matrix  determining  the  coefficients  of  the  polynomial  with  argument  x  formed 
%  from  the  polynomial  PM  of  degree  M  having  as  argument  the  first  degree  polynomial 
%  PI.  Those  coefficients  are  (Pl(2).A(0:M))*A.*(Pl(l).Afliplr(0:M)); 

M^length(PM)-l ; 

A=(ones(M+l,l)*[PM  0])’; 

A~triu(fliplr(pascal(M+ 1  ))).*reshape(  A(  1  :(M  + 1  )A2),M+ 1 JM+ 1 )’; 
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